Introduction

This vignette contains the code used for the case studies described in ‘Incorporating unobserved heterogeneity in Weibull survival models: A Bayesian approach’. All analyses were performed in R using the RMWreg library, which accompanies this manuscript.

Before starting the analysis, we proceed to open a new session of R Studio and to clean the existing R environment using

rm(list = ls())

As a second step, we install and load the RMWreg.

library(RMWreg)
library(parallel) # For parallel run of multiple MCMC chains
library(coda) # To assess convergence of MCMC chains
library(gplots) # To us 'plotCI'
library(compiler); enableJIT(3) # For faster loops
## [1] 0
library(data.table)

In the meantime, the RMWreg can be manually installed using the provided source file. However, the installation process will be simplified once this library is available through the CRAN repository.

Analysis of Autologous and Allogeneic Bone Marrow Transplant dataset

This dataset is publicly available as part of the KMsurv library. It contains post-surgery disease-free survival times (until relapse or death, in months) for 101 advanced acute myelogenous leukemia patients (51 right-censored observations). In the trial, 51 patients received an autologous bone marrow transplant, replacing the patient’s marrow with their own marrow treated with high doses of chemotherapy.The remaining patients had an allogeneic bone transplant, using marrow extracted from a sibling. Only the type of treatment is available as a covariate and thus an important amount of unobserved heterogeneity is expected.

Loading the data and creation of design matrix

library(KMsurv)
data(alloauto)

# Survival time and censoring indicator
Time = alloauto$time # Survival time
Event = alloauto$delta # 1: if observed event; 0: if censored

# Design matrix construction
x0 = rep(1, times = nrow(alloauto)) # Intercept
x1 = alloauto$type-1 # Treatment type - 0: allogeneic; 1: autologous
DesignMat = cbind(x0, x1) 

These data contain subjects.

Descriptive analysis

Number of subjects according to treatment type and censoring status

table(Event, x1)
##      x1
## Event  0  1
##     0 28 23
##     1 22 28

Distribution of follow-up times (relapse, death or censoring) by treatment type

boxplot(Time ~ x1, 
        names = c("Allogeneic", "Autologous"),
        xlab = "Treatment type", ylab = "Follow-up time",
        frame = FALSE)

Distribution of observed event times (relapse or death, excluding censored observations) by treatment type

boxplot(Time[Event == 1] ~ x1[Event == 1], 
        names = c("Allogeneic", "Autologous"),
        xlab = "Treatment type", ylab = "Observed event time",
        frame = FALSE)

Fitting the model

Here, we analyse the data described above using Weibull model (no mixing) and RME-AFT regression models (i.e. RMW-AFT with \(\gamma = 0\)) using all the mixing distributions described in Table 1 of the main manuscript. The following function is used in order to run multiple chains in parallel.

RunRMWreg_MCMC <- function(Row, ParamTable, 
                           Mixing, BaseModel,
                           N, Thin, Burn, Time, Event, DesignMat, 
                           PrintProgress, StoreAdapt, lambdaPeriod,
                           StoreChains = FALSE, StoreDir = getwd(),
                           DataName)
{
  require(RMWreg)
  
  set.seed(ParamTable$seed[Row])
  Chain <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat, 
                       Mixing = Mixing, BaseModel = BaseModel,
                       PrintProgress = PrintProgress,
                       PriorCV = as.character(ParamTable$PriorCV[Row]), 
                       PriorMeanCV = ParamTable$PriorMeanCV[Row],
                       StoreAdapt = StoreAdapt,
                       lambdaPeriod = lambdaPeriod,
                       StoreChains = StoreChains, 
                       RunName = paste0(DataName,"_",Mixing,Row),
                       StoreDir = StoreDir)
  return(Chain)
  
}

Additionally, the following function is used to run graphical and formal convergence diagnostics for MCMC chains

PlotDiag <- function(Chain, titles)
{ 
  for(i in 1:ncol(Chain))
  {
    plot(Chain[,i], main = titles[i], 
         cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5)
  }
  cat("Geweke diagnostics Z scores \n")
  print(round(geweke.diag(Chain)$z,2)); cat("\n") 
  cat("Heidelberger and Welch's convergence diagnostic \n")
  print(round(heidel.diag(Chain)[,2:3],2)); cat("\n") 
  cat("Effective Sample Size \n")
  print(round(effectiveSize(Chain))); cat("\n") 
}

Throughout the following parameters are used

N = 600000; Thin = 50; Burn = 150000
StoreDir = "~/Documents/OneDrive/Projects/Survival/RMW/SuppMaterial/Chains"
titles = c(expression(beta[0]), expression(beta[1]), expression(gamma), expression(theta))

Exponential (no mixing)

ParamTableNone = data.frame(seed = 1)
set.seed(ParamTableNone$seed)
ChainNone <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat, 
                         Mixing = "None", BaseModel = "Exponential",
                         PrintProgress = FALSE, StoreAdapt = TRUE,
                         StoreChains = TRUE, StoreDir = StoreDir,
                         RunName = "AA_None")
ChainNone_mcmc = mcmc(ChainNone$beta)
colnames(ChainNone_mcmc) = titles[1:2]
PlotDiag(ChainNone_mcmc, titles = titles[1:2])

## Geweke diagnostics Z scores 
## beta[0] beta[1] 
##    0.26   -0.45 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.85
## beta[1]     1   0.65
## 
## Effective Sample Size 
## beta[0] beta[1] 
##    9000    9000

RME Exponential mixing

ParamTableExp = data.frame(seed = 2)
set.seed(ParamTableExp$seed)
ChainExp <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat, 
                        Mixing = "Exponential", BaseModel = "Exponential",
                        PrintProgress = FALSE, StoreAdapt = TRUE,
                        StoreChains = TRUE, StoreDir = StoreDir,
                        RunName = "AA_Exp", lambdaPeriod = 1)
ChainExp_mcmc = mcmc(ChainExp$beta)
colnames(ChainExp_mcmc) = titles[1:2]
PlotDiag(ChainExp_mcmc, titles = titles[1:2])

## Geweke diagnostics Z scores 
## beta[0] beta[1] 
##   -0.51   -0.25 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.37
## beta[1]     1   0.07
## 
## Effective Sample Size 
## beta[0] beta[1] 
##    8062    9268

RME Gamma mixing

ParamTableGam <- expand.grid("PriorCV" = c("TruncExp", "Pareto"), 
                             "PriorMeanCV" = c(1.25, 1.5, 2, 5, 10))
set.seed(3)
ParamTableGam$seed = sample(1:1e6, size = nrow(ParamTableGam))
cat("Parameter values \n")
## Parameter values
print(ParamTableGam); cat("\n")
##     PriorCV PriorMeanCV   seed
## 1  TruncExp        1.25 168042
## 2    Pareto        1.25 807516
## 3  TruncExp        1.50 384942
## 4    Pareto        1.50 327734
## 5  TruncExp        2.00 602099
## 6    Pareto        2.00 604392
## 7  TruncExp        5.00 124633
## 8    Pareto        5.00 294599
## 9  TruncExp       10.00 577606
## 10   Pareto       10.00 630974
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainGam <- parLapply(cl, as.list(1:nrow(ParamTableGam)), 
                      fun = RunRMWreg_MCMC,
                      ParamTable = ParamTableGam, 
                      Mixing = "Gamma", BaseModel = "Exponential",
                      N = N, Thin = Thin, Burn = Burn, 
                      Time = Time, Event = Event, DesignMat = DesignMat, 
                      PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
                      StoreChains = TRUE, StoreDir = StoreDir, DataName = "AA")
stopCluster(cl)
ChainGam_mcmc = lapply(ChainGam, function(x) {mcmc(cbind(x$beta, x$theta))})
ChainGam_mcmc = lapply(ChainGam_mcmc, function(x)
                                      {
                                        colnames(x) <- titles[c(1:2,4)]
                                        return(x)
                                      })
Plot = lapply(ChainGam_mcmc, PlotDiag, titles = titles[c(1:2,4)])

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    2.54   -1.64    0.57 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.08
## beta[1]     1   0.72
## theta       1   0.94
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    7296    8346     174

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -1.70    0.22   -0.88 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.30
## beta[1]     1   0.96
## theta       1   0.84
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    5625    9000      41

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.28   -0.77   -0.74 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.59
## beta[1]     1   0.62
## theta       1   0.90
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    5184    9000      39

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.27    0.22   -1.03 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.36
## beta[1]     1   0.59
## theta       1   0.82
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    7633    9000     918

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.15    0.78   -0.84 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.95
## beta[1]     1   0.29
## theta       1   0.71
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    5558    9000     115

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.53    0.79    1.34 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.39
## beta[1]     1   0.63
## theta    2701   0.42
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    8126    9000    1072

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.88    0.05   -1.20 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.08
## beta[1]     1   0.33
## theta       1   0.08
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    8643    8705    1550

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.19    0.35   -0.63 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.80
## beta[1]     1   0.51
## theta       1   0.99
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    8757    9000    2180

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    1.17   -0.42    0.53 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]  3601   0.11
## beta[1]     1   0.49
## theta       1   0.13
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    4535    7861      89

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    0.14    0.52   -0.39 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.95
## beta[1]     1   0.38
## theta       1   0.82
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    7950    9000    1074

RME Inverse Gamma mixing

ParamTableIGam <- expand.grid("PriorCV" = c("TruncExp", "Pareto"), 
                              "PriorMeanCV" = c(1.25, 1.5))
set.seed(4)
ParamTableIGam$seed = sample(1:1e6, size = nrow(ParamTableIGam))
cat("Parameter values \n")
## Parameter values
print(ParamTableIGam); cat("\n")
##    PriorCV PriorMeanCV   seed
## 1 TruncExp        1.25 585801
## 2   Pareto        1.25   8946
## 3 TruncExp        1.50 293740
## 4   Pareto        1.50 277375
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainIGam <- parLapply(cl, as.list(1:nrow(ParamTableIGam)), 
                      fun = RunRMWreg_MCMC,
                      ParamTable = ParamTableIGam, 
                      Mixing = "InvGamma", BaseModel = "Exponential",
                      N = N, Thin = Thin, Burn = Burn, 
                      Time = Time, Event = Event, DesignMat = DesignMat, 
                      PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 5,
                      StoreChains = TRUE, StoreDir = StoreDir, DataName = "AA")
stopCluster(cl)
ChainIGam_mcmc = lapply(ChainIGam, function(x) {mcmc(cbind(x$beta, x$theta))})
ChainIGam_mcmc = lapply(ChainIGam_mcmc, function(x)
                                       {
                                         colnames(x) <- titles[c(1:2,4)]
                                         return(x)
                                       })
Plot = lapply(ChainIGam_mcmc, PlotDiag, titles = titles[c(1:2,4)])

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    1.34   -1.02   -1.32 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.38
## beta[1]     1   0.34
## theta       1   0.53
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    1207    8238     609

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.42   -0.84   -0.29 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.80
## beta[1]     1   0.96
## theta       1   0.94
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##     600    8264     180

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -1.42    0.19    0.99 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.26
## beta[1]     1   0.77
## theta     901   0.68
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    1371    8254     463

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    1.03   -0.38   -1.10 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.38
## beta[1]     1   0.82
## theta       1   0.22
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    1245    8592     418

RME Inverse Gaussian mixing

ParamTableIGauss <- expand.grid("PriorCV" = c("TruncExp", "Pareto"), 
                                "PriorMeanCV" = c(1.25, 1.5, 2))
set.seed(8)
ParamTableIGauss$seed = sample(1:1e6, size = nrow(ParamTableIGauss))
cat("Parameter values \n")
## Parameter values
print(ParamTableIGauss); cat("\n")
##    PriorCV PriorMeanCV   seed
## 1 TruncExp        1.25 466296
## 2   Pareto        1.25 207824
## 3 TruncExp        1.50 799657
## 4   Pareto        1.50 651870
## 5 TruncExp        2.00 321508
## 6   Pareto        2.00 718924
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainIGauss <- parLapply(cl, as.list(1:nrow(ParamTableIGauss)), 
                      fun = RunRMWreg_MCMC,
                      ParamTable = ParamTableIGauss, 
                      Mixing = "InvGauss", BaseModel = "Exponential",
                      N = N, Thin = Thin, Burn = Burn, 
                      Time = Time, Event = Event, DesignMat = DesignMat, 
                      PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 5,
                      StoreChains = TRUE, StoreDir = StoreDir, DataName = "AA")
stopCluster(cl)
ChainIGauss_mcmc = lapply(ChainIGauss, function(x) {mcmc(cbind(x$beta, x$theta))})
ChainIGauss_mcmc = lapply(ChainIGauss_mcmc, function(x)
                                       {
                                         colnames(x) <-  titles[c(1:2,4)]
                                         return(x)
                                       })
Plot = lapply(ChainIGauss_mcmc, PlotDiag, titles = titles[c(1:2,4)])

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.57    0.56   -0.73 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.69
## beta[1]     1   0.50
## theta       1   0.45
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##     870    7752    1150

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    0.12   -0.60   -1.02 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.86
## beta[1]     1   0.66
## theta       1   0.91
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##     448    7242      64

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    0.22   -0.38   -0.99 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.31
## beta[1]     1   0.99
## theta       1   0.14
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    1356    7032     391

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.58   -0.51   -1.24 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]  3601   0.36
## beta[1]     1   0.80
## theta       1   0.05
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##     394    7570    1273

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    0.61    0.74    0.65 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.84
## beta[1]     1   0.74
## theta       1   0.92
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    2062    7431    1886

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    0.94   -2.17    0.59 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.36
## beta[1]     1   0.30
## theta       1   0.64
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    1709    7579    1327

RME Log-Normal mixing

ParamTableLN <- expand.grid("PriorCV" = c("TruncExp", "Pareto"), 
                            "PriorMeanCV" = c(1.25, 1.5, 2, 5, 10))
set.seed(6)
ParamTableLN$seed = sample(1:1e6, size = nrow(ParamTableLN))
cat("Parameter values \n")
## Parameter values
print(ParamTableLN); cat("\n")
##     PriorCV PriorMeanCV   seed
## 1  TruncExp        1.25 606269
## 2    Pareto        1.25 937642
## 3  TruncExp        1.50 264352
## 4    Pareto        1.50 380093
## 5  TruncExp        2.00 807481
## 6    Pareto        2.00 978071
## 7  TruncExp        5.00 957928
## 8    Pareto        5.00 762727
## 9  TruncExp       10.00 509645
## 10   Pareto       10.00  64477
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainLN <- parLapply(cl, as.list(1:nrow(ParamTableLN)), 
                     fun = RunRMWreg_MCMC,
                     ParamTable = ParamTableLN, 
                     Mixing = "LogNormal", BaseModel = "Exponential",
                     N = N, Thin = Thin, Burn = Burn, 
                     Time = Time, Event = Event, DesignMat = DesignMat, 
                     PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 5,
                     StoreChains = TRUE, StoreDir = StoreDir, DataName = "AA")
stopCluster(cl)
ChainLN_mcmc = lapply(ChainLN, function(x) {mcmc(cbind(x$beta, x$theta))})
ChainLN_mcmc = lapply(ChainLN_mcmc, function(x)
                                       {
                                         colnames(x) <- titles[c(1:2,4)]
                                         return(x)
                                       })
Plot = lapply(ChainLN_mcmc, PlotDiag, titles = titles[c(1:2,4)])

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.44    1.41    1.36 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.81
## beta[1]     1   0.09
## theta      NA   0.02
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    7647    6243     946

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    1.62   -2.04   -0.58 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.54
## beta[1]   901   0.11
## theta       1   0.17
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    6217    5807    1169

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    0.12   -0.04   -0.94 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.78
## beta[1]     1   0.62
## theta       1   0.58
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    6860    6638    2093

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.12   -0.48   -0.02 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.47
## beta[1]     1   0.74
## theta       1   0.83
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    5942    5994    1463

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.91    1.36   -0.22 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.80
## beta[1]     1   0.78
## theta       1   0.27
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    5570    5478    2284

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.32    0.73    0.05 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.69
## beta[1]     1   0.13
## theta       1   0.55
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    5105    5001    1548

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.91    0.32    0.15 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.78
## beta[1]     1   0.59
## theta       1   0.14
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    4923    4847    2589

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##    0.99   -1.39    0.82 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.23
## beta[1]     1   0.74
## theta       1   0.53
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    4741    4810    1893

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.26    0.92   -0.56 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.29
## beta[1]     1   0.12
## theta       1   0.67
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    4317    4617    2802

## Geweke diagnostics Z scores 
## beta[0] beta[1]   theta 
##   -0.77    2.19    0.64 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.34
## beta[1]     1   0.29
## theta       1   0.16
## 
## Effective Sample Size 
## beta[0] beta[1]   theta 
##    4339    5005    1659

Weibull (no mixing)

ParamTableWei = data.frame(seed = 7)
set.seed(ParamTableWei$seed)
ChainWei <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat, 
                        Mixing = "None", BaseModel = "Weibull",
                        Hyp1Gam = 4, Hyp2Gam = 1,
                        PrintProgress = FALSE, StoreAdapt = TRUE,
                        StoreChains = TRUE, StoreDir = StoreDir,
                        RunName = "AA_Wei")
ChainWei_mcmc = mcmc(cbind(ChainWei$beta, ChainWei$gam))
colnames(ChainWei_mcmc) = titles[1:3]
PlotDiag(ChainWei_mcmc, titles = titles[1:3])

## Geweke diagnostics Z scores 
## beta[0] beta[1]   gamma 
##   -0.68    1.08    0.16 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.85
## beta[1]     1   0.92
## gamma       1   0.90
## 
## Effective Sample Size 
## beta[0] beta[1]   gamma 
##    9000    8660    9368

RMW Exponential mixing

ParamTableRMWexp = data.frame(seed = 8)
set.seed(ParamTableRMWexp$seed)
ChainRMWexp <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat, 
                           Mixing = "Exponential", BaseModel = "Weibull",
                           Hyp1Gam = 4, Hyp2Gam = 1,
                           PrintProgress = FALSE, StoreAdapt = TRUE,
                           StoreChains = TRUE, StoreDir = StoreDir,
                           RunName = "AA_RMWexp", lambdaPeriod = 1)
ChainRMWexp_mcmc = mcmc(cbind(ChainRMWexp$beta, ChainRMWexp$gam))
colnames(ChainRMWexp_mcmc) = titles[1:3]
PlotDiag(ChainRMWexp_mcmc, titles[1:3])

## Geweke diagnostics Z scores 
## beta[0] beta[1]   gamma 
##    0.23    0.29   -0.23 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.31
## beta[1]     1   0.07
## gamma       1   0.77
## 
## Effective Sample Size 
## beta[0] beta[1]   gamma 
##    8472    8379   10191

Model comparison

Pseudo ML

CDnone = RMWreg_CaseDeletion(ChainNone, Time, Event, DesignMat, 
                             Mixing = "None", BaseModel = "Exponential")
ParamTableNone$logPsML <- sum(CDnone[,1])

CDexp = RMWreg_CaseDeletion(ChainExp, Time, Event, DesignMat, 
                            Mixing = "Exponential", BaseModel = "Exponential")
ParamTableExp$logPsML <- sum(CDexp[,1])

CDgam = lapply(ChainGam, FUN = RMWreg_CaseDeletion, 
               Time = Time, Event = Event, DesignMat = DesignMat, 
               Mixing = "Gamma", BaseModel = "Exponential")
ParamTableGam$logPsML = sapply(CDgam, FUN = function(x) sum(x[,1]))

CDigam = lapply(ChainIGam, FUN = RMWreg_CaseDeletion, 
                Time = Time, Event = Event, DesignMat = DesignMat, 
                Mixing = "InvGamma", BaseModel = "Exponential")
ParamTableIGam$logPsML = sapply(CDigam, FUN = function(x) sum(x[,1]))

CDigauss = lapply(ChainIGauss, FUN = RMWreg_CaseDeletion, 
                  Time = Time, Event = Event, DesignMat = DesignMat, 
                  Mixing = "InvGauss", BaseModel = "Exponential")
ParamTableIGauss$logPsML = sapply(CDigauss, FUN = function(x) sum(x[,1]))

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
CDln = parLapply(cl, ChainLN, fun = RMWreg_CaseDeletion, 
                 Time = Time, Event = Event, DesignMat = DesignMat, 
                 Mixing = "LogNormal", BaseModel = "Exponential")
stopCluster(cl)
ParamTableLN$logPsML = sapply(CDln, FUN = function(x) sum(x[,1]))

CDwei = RMWreg_CaseDeletion(ChainWei, Time, Event, DesignMat, 
                            Mixing = "None", BaseModel = "Weibull")
ParamTableWei$logPsML <- sum(CDwei[,1])

CDrmwexp = RMWreg_CaseDeletion(ChainRMWexp, Time, Event, DesignMat, 
                               Mixing = "Exponential", BaseModel = "Weibull")
ParamTableRMWexp$logPsML <- sum(CDrmwexp[,1])

ML

RunRMWreg_logML <- function(Row, ParamTable, Chains,
                            Mixing, BaseModel,
                            Time, Event, DesignMat, 
                            Hyp1Gam = 1, Hyp2Gam = 1, N, Thin, Burn)
{
  require(RMWreg)
  
  logML = RMWreg_logML(Chains[[Row]], Time, Event, DesignMat, 
                       PriorCV = as.character(ParamTable$PriorCV[Row]),
                       PriorMeanCV = ParamTable$PriorMeanCV[Row],
                       Hyp1Gam = Hyp1Gam, Hyp2Gam = Hyp2Gam,
                       Mixing = Mixing, BaseModel = BaseModel, N = N, Thin = Thin, Burn = Burn)
  
  return(logML)
  
}

ParamTableNone$logML = RMWreg_logML(ChainNone, Time, Event, DesignMat, 
                                    Mixing = "None", BaseModel = "Exponential", 
                                    N = N, Thin = Thin, Burn = Burn)

ParamTableExp$logML = RMWreg_logML(ChainExp, Time, Event, DesignMat, 
                                   Mixing = "Exponential", BaseModel = "Exponential", 
                                   N = N, Thin = Thin, Burn = Burn)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableGam$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableGam)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableGam, 
                                      Chains = ChainGam,
                                      Mixing = "Gamma", BaseModel = "Exponential",
                                      Time = Time, Event = Event, DesignMat = DesignMat, 
                                      N = N, Thin = Thin, Burn = Burn))
stopCluster(cl)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableIGam$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableIGam)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableIGam, 
                                      Chains = ChainIGam,
                                      Mixing = "InvGamma", BaseModel = "Exponential",
                                      Time = Time, Event = Event, DesignMat = DesignMat, 
                                      N = N, Thin = Thin, Burn = Burn))
stopCluster(cl)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableIGauss$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableIGauss)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableIGauss, 
                                      Chains = ChainIGauss,
                                      Mixing = "InvGauss", BaseModel = "Exponential",
                                      Time = Time, Event = Event, DesignMat = DesignMat, 
                                      N = N, Thin = Thin, Burn = Burn))
stopCluster(cl)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableLN$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableLN)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableLN, 
                                      Chains = ChainLN,
                                      Mixing = "LogNormal", BaseModel = "Exponential",
                                      Time = Time, Event = Event, DesignMat = DesignMat, 
                                      N = N, Thin = Thin, Burn = Burn))
stopCluster(cl)

ParamTableWei$logML = RMWreg_logML(ChainWei, Time, Event, DesignMat, 
                                   Mixing = "None", BaseModel = "Weibull",
                                   Hyp1Gam = 4, Hyp2Gam = 1, 
                                   N = N, Thin = Thin, Burn = Burn)

ParamTableRMWexp$logML = RMWreg_logML(ChainWei, Time, Event, DesignMat, 
                                      Mixing = "Exponential", BaseModel = "Weibull",
                                      Hyp1Gam = 4, Hyp2Gam = 1, 
                                      N = N, Thin = Thin, Burn = Burn)

DIC

ParamTableNone$DIC <- RMWreg_DIC(ChainNone, Time, Event, DesignMat, 
                                 Mixing = "None", BaseModel = "Exponential")

ParamTableExp$DIC <- RMWreg_DIC(ChainExp, Time, Event, DesignMat, 
                                Mixing = "Exponential", BaseModel = "Exponential")

ParamTableGam$DIC = sapply(ChainGam, FUN = RMWreg_DIC, 
                           Time = Time, Event = Event, DesignMat = DesignMat, 
                           Mixing = "Gamma", BaseModel = "Exponential")

ParamTableIGam$DIC = sapply(ChainIGam, FUN = RMWreg_DIC, 
                            Time = Time, Event = Event, DesignMat = DesignMat, 
                            Mixing = "InvGamma", BaseModel = "Exponential")

ParamTableIGauss$DIC = sapply(ChainIGauss, FUN = RMWreg_DIC, 
                              Time = Time, Event = Event, DesignMat = DesignMat, 
                              Mixing = "InvGauss", BaseModel = "Exponential")

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableLN$DIC = simplify2array(parLapply(cl, ChainLN, fun = RMWreg_DIC, 
                                            Time = Time, Event = Event, DesignMat = DesignMat, 
                                            Mixing = "LogNormal", BaseModel = "Exponential"))
stopCluster(cl)

ParamTableWei$DIC <- RMWreg_DIC(ChainWei, Time, Event, DesignMat, 
                                Mixing = "None", BaseModel = "Weibull")

ParamTableRMWexp$DIC <- RMWreg_DIC(ChainWei, Time, Event, DesignMat, 
                                   Mixing = "Exponential", BaseModel = "Weibull")

Summary

cat("\n RME no mixing \n")
## 
##  RME no mixing
ParamTableNone
##   seed   logPsML     logML      DIC
## 1    1 -230.4599 -229.4623 459.9963
cat("\n RME Exponential mixing \n")
## 
##  RME Exponential mixing
ParamTableExp
##   seed   logPsML     logML      DIC
## 1    2 -223.5881 -222.0672 446.5538
cat("\n RME Gamma mixing \n")
## 
##  RME Gamma mixing
ParamTableGam
##     PriorCV PriorMeanCV   seed   logPsML     logML      DIC
## 1  TruncExp        1.25 168042 -228.3560 -227.4952 455.9293
## 2    Pareto        1.25 807516 -228.2476 -227.0317 455.7122
## 3  TruncExp        1.50 384942 -227.6544 -226.7201 454.6764
## 4    Pareto        1.50 327734 -227.5672 -226.8410 454.4988
## 5  TruncExp        2.00 602099 -227.2252 -226.1494 453.8678
## 6    Pareto        2.00 604392 -227.2272 -226.4229 453.8926
## 7  TruncExp        5.00 124633 -226.6552 -225.3898 452.6433
## 8    Pareto        5.00 294599 -226.9003 -225.9342 453.2963
## 9  TruncExp       10.00 577606 -226.5228 -225.2106 452.3818
## 10   Pareto       10.00 630974 -226.8806 -225.8763 453.2428
cat("\n RME Inverse Gamma mixing \n")
## 
##  RME Inverse Gamma mixing
ParamTableIGam
##    PriorCV PriorMeanCV   seed   logPsML     logML      DIC
## 1 TruncExp        1.25 585801 -224.5198 -224.6049 448.8781
## 2   Pareto        1.25   8946 -224.5838 -224.7460 449.0500
## 3 TruncExp        1.50 293740 -224.2709 -223.8656 448.3744
## 4   Pareto        1.50 277375 -224.2931 -224.0833 448.4058
cat("\n RME Inverse Gaussian mixing \n")
## 
##  RME Inverse Gaussian mixing
ParamTableIGauss
##    PriorCV PriorMeanCV   seed   logPsML     logML      DIC
## 1 TruncExp        1.25 466296 -223.8578 -224.2714 447.7325
## 2   Pareto        1.25 207824 -223.9478 -224.2283 448.0506
## 3 TruncExp        1.50 799657 -223.4771 -223.2344 447.0858
## 4   Pareto        1.50 651870 -223.6137 -223.4765 447.3825
## 5 TruncExp        2.00 321508 -223.3694 -222.6429 446.8220
## 6   Pareto        2.00 718924 -223.3786 -222.9020 446.8793
cat("\n RME Log-Normal mixing \n")
## 
##  RME Log-Normal mixing
ParamTableLN
##     PriorCV PriorMeanCV   seed   logPsML     logML      DIC
## 1  TruncExp        1.25 606269 -225.2078 -225.5257 449.9074
## 2    Pareto        1.25 937642 -224.3059 -225.0513 448.2485
## 3  TruncExp        1.50 264352 -223.7999 -223.7006 447.3188
## 4    Pareto        1.50 380093 -223.4678 -223.5243 446.8311
## 5  TruncExp        2.00 807481 -223.0690 -222.3484 446.0267
## 6    Pareto        2.00 978071 -223.2085 -222.6695 446.3980
## 7  TruncExp        5.00 957928 -222.7849 -221.2237 445.6284
## 8    Pareto        5.00 762727 -223.0131 -222.0602 446.0679
## 9  TruncExp       10.00 509645 -222.9367 -221.2571 445.9855
## 10   Pareto       10.00  64477 -223.0865 -221.9649 446.2496
cat("\n RMW no mixing \n")
## 
##  RMW no mixing
ParamTableWei
##   seed   logPsML    logML      DIC
## 1    7 -225.0809 -227.927 450.1554
cat("\n RMW Exponential mixing \n")
## 
##  RMW Exponential mixing
ParamTableRMWexp
##   seed   logPsML     logML    DIC
## 1    8 -223.3197 -224.6084 450.69

Figure 5

# Comparison Weibull vs Exponential
round(exp(ParamTableWei$logML - ParamTableNone$logML), 1)
## [1] 4.6
# Posterior summary for gam under Weibull model
round(median(ChainWei$gam),2)
## [1] 0.69
round(HPDinterval(mcmc(ChainWei$gam)),2)
##    lower upper
## V1  0.54  0.87
## attr(,"Probability")
## [1] 0.95
# Comparison RMWexp vs RMEexp
round(exp(ParamTableExp$logML - ParamTableRMWexp$logML), 1)
## [1] 12.7
# Posterior summary for gam under Weibull model
round(median(ChainRMWexp$gam),2)
## [1] 0.86
round(HPDinterval(mcmc(ChainRMWexp$gam)),2)
##    lower upper
## V1  0.65  1.06
## attr(,"Probability")
## [1] 0.95
# log-BFs
LOGBF1 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
           ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 1.25],
           ParamTableIGam$logML[ParamTableIGam$PriorCV == "TruncExp" & ParamTableIGam$PriorMeanCV == 1.25],
           ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 1.25],
           ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 1.25]) -
          ParamTableNone$logML
LOGBF1_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 1.25],
             ParamTableIGam$logML[ParamTableIGam$PriorCV == "Pareto" & ParamTableIGam$PriorMeanCV == 1.25],
             ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 1.25],
             ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 1.25]) - 
           ParamTableNone$logML
LOGBF2 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
           ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 1.5],
           ParamTableIGam$logML[ParamTableIGam$PriorCV == "TruncExp" & ParamTableIGam$PriorMeanCV == 1.5],
           ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 1.5],
           ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 1.5]) -
          ParamTableNone$logML
LOGBF2_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 1.5],
             ParamTableIGam$logML[ParamTableIGam$PriorCV == "Pareto" & ParamTableIGam$PriorMeanCV == 1.5],
             ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 1.5],
             ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 1.5]) - 
           ParamTableNone$logML

LOGBF3 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
           ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 2],
           ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 2],
           ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 2]) -
          ParamTableNone$logML
LOGBF3_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 2],
             ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 2],
             ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 2]) - 
           ParamTableNone$logML

LOGBF4 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
           ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 5],
           ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 5]) -
          ParamTableNone$logML
LOGBF4_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 5],
             ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 5]) - 
           ParamTableNone$logML

LOGBF5 = c(ParamTableNone$logML, ParamTableWei$logML, ParamTableExp$logML,
           ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 10],
           ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 10]) -
          ParamTableNone$logML
LOGBF5_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 10],
             ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 10]) - 
           ParamTableNone$logML

# log-PsBFs
LOGPsBF1 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
           ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 1.25],
           ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "TruncExp" & ParamTableIGam$PriorMeanCV == 1.25],
           ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 1.25],
           ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 1.25]) -
          ParamTableNone$logPsML
LOGPsBF1_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 1.25],
             ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "Pareto" & ParamTableIGam$PriorMeanCV == 1.25],
             ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 1.25],
             ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 1.25]) - 
           ParamTableNone$logPsML
LOGPsBF2 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
           ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 1.5],
           ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "TruncExp" & ParamTableIGam$PriorMeanCV == 1.5],
           ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 1.5],
           ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 1.5]) -
          ParamTableNone$logPsML
LOGPsBF2_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 1.5],
             ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "Pareto" & ParamTableIGam$PriorMeanCV == 1.5],
             ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 1.5],
             ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 1.5]) - 
           ParamTableNone$logPsML

LOGPsBF3 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
           ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 2],
           ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "TruncExp" & ParamTableIGauss$PriorMeanCV == 2],
           ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 2]) -
          ParamTableNone$logPsML
LOGPsBF3_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 2],
             ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "Pareto" & ParamTableIGauss$PriorMeanCV == 2],
             ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 2]) - 
           ParamTableNone$logPsML

LOGPsBF4 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
           ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 5],
           ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 5]) -
          ParamTableNone$logPsML
LOGPsBF4_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 5],
             ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 5]) - 
           ParamTableNone$logPsML

LOGPsBF5 = c(ParamTableNone$logPsML, ParamTableWei$logPsML, ParamTableExp$logPsML,
           ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" & ParamTableGam$PriorMeanCV == 10],
           ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & ParamTableLN$PriorMeanCV == 10]) -
          ParamTableNone$logPsML
LOGPsBF5_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" & ParamTableGam$PriorMeanCV == 10],
             ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & ParamTableLN$PriorMeanCV == 10]) - 
           ParamTableNone$logPsML


par(mfrow=c(2,3))
par(mar=c(5,5,4,0.5))
plot(LOGBF1, LOGPsBF1, ylim = c(0,8.5), xlim = c(0,8.5),
     ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
     pch = c(4,3,8,0,5,1,2), cex = rep(3,7), bty = "n",
     main = expression(paste("E(",c[v],")=1.25")), 
     cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
points(LOGBF1_1, LOGPsBF1_1, pch = c(15,18,16,17), lwd = 2, cex = c(3,4.5,3,3))
plot(LOGBF2, LOGPsBF2, ylim = c(0,8.5), xlim = c(0,8.5),
     ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
     pch = c(4,3,8,0,5,1,2), cex = rep(3,7), bty = "n",
     main = expression(paste("E(",c[v],")=1.5")), 
     cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd =2)
points(LOGBF2_1, LOGPsBF2_1, pch = c(15,18,16,17), lwd = 2, cex = c(3,4.5,3,3))
plot(LOGBF3, LOGPsBF3, ylim = c(0,8.5), xlim = c(0,8.5), 
     ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
     pch = c(4,3,8,0,1,2), cex = rep(3,6), bty = "n",
     main = expression(paste("E(",c[v],")=2")), 
     cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
points(LOGBF3_1, LOGPsBF3_1, pch = c(15,16,17), lwd = 2, cex = rep(3,4))
plot(LOGBF4, LOGPsBF4, ylim = c(0,8.5), xlim = c(0,8.5),
     ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
     pch = c(4,3,8,0,2), cex = rep(3,5), bty = "n",
     main = expression(paste("E(",c[v],")=5")),
     cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
points(LOGBF4_1, LOGPsBF4_1, pch = c(15,17), lwd = 2, cex = rep(3,4))
plot(LOGBF5, LOGPsBF5, ylim = c(0,8.5), xlim = c(0,8.5),
     ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors",
     pch = c(4,3,8,0,2), cex = rep(3,5), bty="n",
     main = expression(paste("E(",c[v],")=10")),
     cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
points(LOGBF5_1, LOGPsBF5_1, pch = c(15,17), lwd = 2, cex = rep(3,4))
plot(1, type = "n", axes = F, xlab = "", ylab = "")
legend('center', c("Exponential","Weibull","RME-exp","RME-gam",
                   "RME-inv. gam","RME-inv. Gauss","RME-log normal"),
       pch = c(4,3,8,0,5,1,2), col = "black", cex = 3, bty = "n")

Figure 6

prior = function(x, a, trunc, upper=Inf, type)
{
    if(type=="TruncExp")
    {
        if(trunc==1){aux=dexp(x,rate=a)*exp(a)}
        if(trunc==2){aux=dexp(x,rate=a)/(exp(-a)-exp(-upper*a))}
    }
    if(type=="Pareto")
    {
        if(trunc==1){aux=a*x^(-(a+1))}
        if(trunc==2){aux=(a*x^(-(a+1)))/(1-upper^(-a))}
    }
    return(aux)
}

# CV as a function of theta
CV.RMEGAM<-function(theta) { sqrt(theta/(theta-2)) }
CV.RMEIGAM<-function(theta) { sqrt((theta+2)/(theta)) }
CV.RMEIGAUSS<-function(theta){ sqrt((5*theta^2+4*theta+1)/(theta^2+2*theta+1)) }
CV.RMELN<-function(theta){ sqrt(2*exp(theta)-1) }

par(mfrow=c(2,2))
par(mar=c(5, 5, 4, 2))

cv = seq(1, 50, by = 0.01)
hist(CV.RMELN(ChainLN[[7]]$theta), freq = FALSE, 
     ylim = c(0,0.4), xlim = c(1,10), breaks = 80, 
     main = expression(paste("Trunc. Exp. prior with ","E(",cv,")=5")),
     ylab = "Frequency/Density", xlab = expression(R[cv]),
     cex.lab = 2, cex.axis = 1.5, cex.main = 2.5)
lines(cv, prior(cv, a = 0.25, trunc = 1, upper = Inf, type = "TruncExp"), lwd = 4)
hist(CV.RMELN(ChainLN[[8]]$theta), freq = FALSE,
     ylim = c(0,0.4), xlim = c(1,10), breaks = 100,
     main = expression(paste("Pareto prior with ","E(",cv,")=5 ")),
     ylab = "Frequency/Density", xlab = expression(R[cv]),
     cex.lab = 2, cex.axis = 1.5, cex.main = 2.5)
lines(cv, prior(cv, a = 1.25, trunc = 1, upper = Inf, type = "Pareto"), lwd = 4)

hist(CV.RMELN(ChainLN[[9]]$theta), freq = FALSE,
     ylim = c(0,0.3), xlim = c(1,10), breaks = 150,
     main = expression(paste("Trunc. Exp. prior with ","E(",cv,")=10")),
     ylab = "Frequency/Density", xlab = expression(R[cv]), 
     cex.lab = 2, cex.axis = 1.5, cex.main = 2.5)
lines(cv, prior(cv, a = 0.10, trunc = 1, upper = Inf, type = "TruncExp"), lwd = 4)

hist(CV.RMELN(ChainLN[[10]]$theta), freq = FALSE,
     ylim = c(0,0.3), xlim = c(1,10), breaks = 200,
     main = expression(paste("Pareto prior with ","E(",cv,")=10")),
     ylab = "Frequency/Density", xlab = expression(R[cv]),
     cex.lab = 2, cex.axis = 1.5, cex.main = 2.5)
lines(cv, prior(cv, a = 1.10, trunc = 1, upper = Inf, type = "Pareto"), lwd = 4)

Figure 7

This figure summarizes posterior estimates for the regression coefficients: posterior median and 95% Highest Posterior Density (HPD) interval

PM.b0 = c(median(ChainNone$beta[,1]), 
          median(ChainWei$beta[,1]),  
          median(ChainExp$beta[,1]),
          unlist(lapply(ChainGam, function(x) median(x$beta[,1]))),
          unlist(lapply(ChainIGam, function(x) median(x$beta[,1]))),
          unlist(lapply(ChainIGauss, function(x) median(x$beta[,1]))),
          unlist(lapply(ChainLN, function(x) median(x$beta[,1]))))
HPD.b0 = rbind(HPDinterval(ChainNone_mcmc[,1]), 
               HPDinterval(ChainWei_mcmc[,1]),
               HPDinterval(ChainExp_mcmc[,1]),
         t(sapply(ChainGam_mcmc, function(x) HPDinterval(x[,1]))),
         t(sapply(ChainIGam_mcmc, function(x) HPDinterval(x[,1]))),
         t(sapply(ChainIGauss_mcmc, function(x) HPDinterval(x[,1]))),
         t(sapply(ChainLN_mcmc, function(x) HPDinterval(x[,1]))))

PM.b1 = c(median(ChainNone$beta[,2]), 
          median(ChainWei$beta[,2]),  
          median(ChainExp$beta[,2]),
          unlist(lapply(ChainGam, function(x) median(x$beta[,2]))),
          unlist(lapply(ChainIGam, function(x) median(x$beta[,2]))),
          unlist(lapply(ChainIGauss, function(x) median(x$beta[,2]))),
          unlist(lapply(ChainLN, function(x) median(x$beta[,2]))))
HPD.b1 = rbind(HPDinterval(ChainNone_mcmc[,2]), 
               HPDinterval(ChainWei_mcmc[,2]), 
               HPDinterval(ChainExp_mcmc[,2]),
               t(sapply(ChainGam_mcmc, function(x) HPDinterval(x[,2]))),
               t(sapply(ChainIGam_mcmc, function(x) HPDinterval(x[,2]))),
               t(sapply(ChainIGauss_mcmc, function(x) HPDinterval(x[,2]))),
               t(sapply(ChainLN_mcmc, function(x) HPDinterval(x[,2]))))

par(mfrow = c(2, 1))
par(mar = c(5,6,4,0.5))

plotCI(x = PM.b0, uiw = HPD.b0[,2] - PM.b0, liw = PM.b0 - HPD.b0[,1],
       main = expression("Intercept"), xaxt = "n", lwd = 2, 
       xlab = "Mixing-Prior", ylab = expression(beta[0]),
       cex.lab = 2, cex.axis = 2, cex.main = 2.5, gap = 1.5)
axis(side = 1, at = c(1.5,3,6,11,14.3,16.6,19,22,26,31),
     labels = c("None", "Exp", "Gamma-Trunc.Exp.", "Gam-Pareto", 
                "IGam-T.E.", "IGam-P", "IGauss-T.E.", "I.Gauss-P.", 
                "LN-Trunc.Exp", "LN-Pareto"), cex.axis = 1)
abline(v = c(2.5,3.5,8.5,13.5,15.5,17.5,20.5,23.5,28.5), lty = 2)
abline(h = 0, lty = 3)

plotCI(x = PM.b1, uiw = HPD.b1[,2] - PM.b1, liw = PM.b1 - HPD.b1[,1],
       main = expression("Treatment: autologous"), xaxt = "n", lwd = 2,
       xlab = "Mixing-Prior", ylab = expression(beta[1]),
       cex.lab = 2, cex.axis = 2, cex.main = 2.5, gap = 1.5)
axis(side = 1, at = c(1.5,3,6,11,14.3,16.6,19,22,26,31),
     labels = c("None", "Exp", "Gamma-Trunc.Exp.", "Gam-Pareto", 
                "IGam-T.E.", "IGam-P", "IGauss-T.E.", "I.Gauss-P.",
                "LN-Trunc.Exp", "LN-Pareto"), cex.axis = 1)
abline(v = c(2.5,3.5,8.5,13.5,15.5,17.5,20.5,23.5,28.5), lty = 2)
abline(h = 0, lty = 3)

Figure 8

BoxplotLambda <- function(Chain, main, ref, ylim, cex.lab, cex.axis, cex.main)
{
    PM = as.vector(apply(Chain$lambda, 2, "median"))
    U = as.vector(HPDinterval(mcmc(Chain$lambda))[,2]) - PM
    L = PM - as.vector(HPDinterval(mcmc(Chain$lambda))[,1])

    plotCI(x = PM, uiw = U, liw = L,
           main = main, xaxt = "n", lwd = 1.5,
           xlab = "Patient", ylab = expression(lambda[i]),
           cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, ylim = ylim, cex = 1)
    abline(h = ref, lty = 2)

    return(PM)
}

# Outlier detection
BFexp = RMWreg_BFoutlier(ChainExp, RefLambda = Event + (1-Event)*1/2,
                         Time, Event, DesignMat,
                         Mixing = "Exponential", BaseModel = "Exponential")

par(mar = c(5,6,4,0.5))
par(mfrow = c(2,1))
PM.RME.EXP = BoxplotLambda(Chain = ChainExp, main = "(a)",
                           ref = c(0.5,1), ylim = c(0,5), cex.lab = 2, cex.axis = 2, cex.main=2)
## Note: no visible binding for global variable 'lambda' 
## Note: no visible binding for global variable 'i'
points((1:length(Time))[Event == 0], PM.RME.EXP[Event == 0], lwd = 2, pch = 16)
axis(side = 1, at = c(25,75), labels = c("Treatment: Allogeneic","Treatment: Autologous"), cex.axis = 2)

plot(BFexp, type = "l", ylim = c(0,2), main = "(b)",
     cex.main = 2, cex.lab = 2, cex.axis = 2, ylab = "Bayes Factors", xlab = "Patient")
abline(h = 1, lty = 2)

Analysis of Cerebral Palsy dataset

These data contain records of 1,549 children affected by cerebral palsy and born during the period 1966-1984 in the area of the Mersey Region Health Authority. The time to follow-up (survival or censoring) is recorded in years since birth. Following the analysis in Kwong and Hutton (2003), we use the amount of severe impairments (ambulation, manual dexterity and mental ability) and the birth weight (in kg.) as predictors for the time to death.

These data were kindly provided by P.O.D. Pharoah and Jane Hutton and we are not authorized to release it. However, this file contains relevant summary information and the code used to perform the analysis. These should serve as a guidance for potential users of the RMWreg library.

Loading the data and creation of design matrix

In this example, the data file has been stored in a data.path folder.

Data = read.csv(file.path(data.path, "CPdata.csv"), sep="")

# Removal of subjects with missing birth weight
Data = Data[is.na(Data$bw)==FALSE,] 

# Survival time and censoring indicator
Time = Data$Lifeyr # Survival time
Event = Data$Dead # 1: if observed event; 0: if censored

# Design matrix construction 
x0 = rep(1, times = nrow(Data)) # Intercept
x1 = Data$Ambul + Data$Mandex + Data$Mentab # No of impairments (treated as categorical)
x1.0 = as.numeric(I(x1 == 0)) 
x1.1 = as.numeric(I(x1 == 1))
x1.2 = as.numeric(I(x1 == 2))
x2 = Data$bw / 1000 # Birth weights in kgs
DesignMat = cbind(x0, x1.0, x1.1, x1.2, x2)

These data contain 1549 subjects.

Descriptive analysis

Number of subjects according to the number of impairments and censoring status

table(Event, x1)
##      x1
## Event   0   1   2   3
##     0 960 208  49  90
##     1  23  25  21 173

Number of subjects according to birth weight and censoring status

bw_group = ifelse(Data$bw < 1500, "1. Very low",
                  ifelse(Data$bw < 2500, "2. Low", "3. Normal"))
table(Event, bw_group)
##      bw_group
## Event 1. Very low 2. Low 3. Normal
##     0         197    341       769
##     1          24     56       162

Distribution of follow-up times (survival or censoring) by the number of impairments

boxplot(Time ~ x1, 
        names = 0:3,
        xlab = "No. of impairments", ylab = "Follow-up time",
        frame = FALSE)

Distribution of observed event times (relapse or death, excluding censored observations) by the number of impairments

boxplot(Time[Event == 1] ~ x1[Event == 1], 
        names = 0:3,
        xlab = "No. of impairments", ylab = "Observed event time",
        frame = FALSE)

Distribution of follow-up times (survival or censoring) by birth weight group

boxplot(Time ~ bw_group, 
        xlab = "Birth weight group", ylab = "Follow-up time",
        frame = FALSE)

Distribution of observed event times (relapse or death, excluding censored observations) by birth weight group

boxplot(Time[Event == 1] ~ bw_group[Event == 1], 
        xlab = "Birth weight group", ylab = "Observed event time",
        frame = FALSE)

RMW-AFT analysis

RunRMWreg_MCMC <- function(Row, ParamTable, 
                           Mixing, BaseModel,
                           N, Thin, Burn, Time, Event, DesignMat, 
                           PrintProgress, StoreAdapt, lambdaPeriod,
                           StoreChains = FALSE, StoreDir = getwd(),
                           DataName)
{
  require(RMWreg)
  
  if("PriorCV" %in% names(ParamTable))
  {
    set.seed(ParamTable$seed[Row])
    Chain <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat, 
                        Mixing = Mixing, BaseModel = BaseModel,
                        PrintProgress = PrintProgress,
                        PriorCV = as.character(ParamTable$PriorCV[Row]), 
                        PriorMeanCV = ParamTable$PriorMeanCV[Row],
                        Hyp1Gam = ParamTable$Hyp1Gam[Row],
                        Hyp2Gam = ParamTable$Hyp2Gam[Row],
                        StoreAdapt = StoreAdapt,
                        lambdaPeriod = lambdaPeriod,
                        StoreChains = StoreChains, 
                        RunName = paste0(DataName,"_",Mixing,Row),
                        StoreDir = StoreDir)
  }
  else
  {
    set.seed(ParamTable$seed[Row])
    Chain <- RMWreg_MCMC(N, Thin, Burn, Time, Event, DesignMat, 
                        Mixing = Mixing, BaseModel = BaseModel,
                        PrintProgress = PrintProgress,
                        Hyp1Gam = ParamTable$Hyp1Gam[Row],
                        Hyp2Gam = ParamTable$Hyp2Gam[Row],
                        StoreAdapt = StoreAdapt,
                        lambdaPeriod = lambdaPeriod,
                        StoreChains = StoreChains, 
                        RunName = paste0(DataName,"_",Mixing,Row),
                        StoreDir = StoreDir)
  }

  return(Chain)
  
}

Additionally, the following function is used to run graphical and formal convergence diagnostics for MCMC chains

PlotDiag <- function(Chain, titles)
{ 
  for(i in 1:ncol(Chain))
  {
    plot(Chain[,i], main = titles[i], 
         cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5)
  }
  cat("Geweke diagnostics Z scores \n")
  print(round(geweke.diag(Chain)$z,2)); cat("\n") 
  cat("Heidelberger and Welch's convergence diagnostic \n")
  print(round(heidel.diag(Chain)[,2:3],2)); cat("\n") 
  cat("Effective Sample Size \n")
  print(round(effectiveSize(Chain))); cat("\n") 
}

Throughout the following parameters are used

N = 600000; Thin = 50; Burn = 150000
StoreDir = "~/Documents/OneDrive/Projects/Survival/RMW/SuppMaterial/Chains"
titles  = c(expression(beta[0]), expression(beta[1]), 
            expression(beta[2]), expression(beta[3]), expression(beta[4]),
            expression(gamma), expression(theta))

Weibull (no mixing)

ParamTableNone <- cbind.data.frame("Hyp1Gam" = c(4, 1), "Hyp2Gam" = c(1, 1))
set.seed(100)
ParamTableNone$seed = sample(1:1e6, size = nrow(ParamTableNone))
cat("Parameter values \n")
## Parameter values
print(ParamTableNone); cat("\n")
##   Hyp1Gam Hyp2Gam   seed
## 1       4       1 307767
## 2       1       1 257673
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainNone <- parLapply(cl, as.list(1:nrow(ParamTableNone)), 
                       fun = RunRMWreg_MCMC,
                       ParamTable = ParamTableNone, 
                       Mixing = "None", BaseModel = "Weibull",
                       N = N, Thin = Thin, Burn = Burn, 
                       Time = Time, Event = Event, DesignMat = DesignMat, 
                       PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
                       StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)
ChainNone_mcmc = lapply(ChainNone, function(x) {mcmc(cbind(x$beta, x$gam))})
ChainNone_mcmc = lapply(ChainNone_mcmc, function(x)
                                      {
                                        colnames(x) <- titles[1:6]
                                        return(x)
                                      })
Plot = lapply(ChainNone_mcmc, PlotDiag, titles = titles[1:6])

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma 
##    1.33    0.19    0.43    1.03   -1.37   -1.01 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.66
## beta[1]     1   0.31
## beta[2]     1   0.94
## beta[3]     1   0.67
## beta[4]     1   0.44
## gamma       1   0.47
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma 
##    3653    9000    9490    8667    3761    8594

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma 
##   -1.06   -0.71    0.74   -0.18    1.02   -0.05 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.55
## beta[1]     1   0.62
## beta[2]     1   0.50
## beta[3]     1   0.52
## beta[4]     1   0.59
## gamma       1   0.26
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma 
##    3975    9000    9000    9000    3997    9000

RMW Exponential mixing

ParamTableExp <- cbind.data.frame("Hyp1Gam" = c(4, 1), "Hyp2Gam" = c(1, 1))
set.seed(101)
ParamTableExp$seed = sample(1:1e6, size = nrow(ParamTableExp))
cat("Parameter values \n")
## Parameter values
print(ParamTableExp); cat("\n")
##   Hyp1Gam Hyp2Gam   seed
## 1       4       1 372199
## 2       1       1  43825
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainExp <- parLapply(cl, as.list(1:nrow(ParamTableExp)), 
                      fun = RunRMWreg_MCMC,
                      ParamTable = ParamTableExp, 
                      Mixing = "Exponential", BaseModel = "Weibull",
                      N = N, Thin = Thin, Burn = Burn, 
                      Time = Time, Event = Event, DesignMat = DesignMat, 
                      PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
                      StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)
ChainExp_mcmc = lapply(ChainExp, function(x) {mcmc(cbind(x$beta, x$gam))})
ChainExp_mcmc = lapply(ChainExp_mcmc, function(x)
                                      {
                                        colnames(x) <- titles[1:6]
                                        return(x)
                                      })
Plot = lapply(ChainExp_mcmc, PlotDiag, titles = titles[1:6])

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma 
##    1.19   -2.67   -0.26    0.25   -0.98    1.17 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.36
## beta[1]   901   0.22
## beta[2]     1   0.76
## beta[3]     1   0.89
## beta[4]     1   0.42
## gamma       1   0.12
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma 
##    2514    9000    9000    9000    2644    9000

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma 
##    0.21    1.55    1.03   -1.82   -0.04   -1.13 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.67
## beta[1]     1   0.49
## beta[2]     1   0.27
## beta[3]   901   0.23
## beta[4]     1   0.61
## gamma       1   0.66
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma 
##    2519    8387    9000    9000    2593    8545

RMW Gamma mixing

ParamTableGam <- expand.grid("Hyp1Gam" = c(4, 1),
                             "PriorMeanCV" = c(1.5, 5),
                             "PriorCV" = c("TruncExp", "Pareto"))
ParamTableGam$Hyp2Gam = with(ParamTableGam, 
                             ifelse(Hyp1Gam == 4, 1, ifelse(Hyp1Gam == 1, 1, 0.1)))
set.seed(102)
ParamTableGam$seed = sample(1:1e6, size = nrow(ParamTableGam))
cat("Parameter values \n")
## Parameter values
print(ParamTableGam); cat("\n")
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed
## 1       4         1.5 TruncExp       1 571629
## 2       1         1.5 TruncExp       1 492883
## 3       4         5.0 TruncExp       1 783694
## 4       1         5.0 TruncExp       1 540340
## 5       4         1.5   Pareto       1  88002
## 6       1         1.5   Pareto       1 420628
## 7       4         5.0   Pareto       1 976328
## 8       1         5.0   Pareto       1 330313
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainGam <- parLapply(cl, as.list(1:nrow(ParamTableGam)), 
                      fun = RunRMWreg_MCMC,
                      ParamTable = ParamTableGam, 
                      Mixing = "Gamma", BaseModel = "Weibull",
                      N = 2*N, Thin = 2*Thin, Burn = 2*Burn, 
                      Time = Time, Event = Event, DesignMat = DesignMat, 
                      PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
                      StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)
ChainGam_mcmc = lapply(ChainGam, function(x) {mcmc(cbind(x$beta, x$gam, x$theta))})
ChainGam_mcmc = lapply(ChainGam_mcmc, function(x)
                                      {
                                        colnames(x) <- titles
                                        return(x)
                                      })
Plot = lapply(ChainGam_mcmc, PlotDiag, titles = titles)

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.04    0.37   -1.09    0.62   -0.12    0.91   -0.41 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.80
## beta[1]     1   0.65
## beta[2]     1   0.69
## beta[3]     1   0.99
## beta[4]     1   0.66
## gamma       1   0.22
## theta       1   0.97
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    5231    7794    8025    9000    5446    4483    1241

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    2.38    2.33    1.00   -0.17   -1.37   -1.71    1.43 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.18
## beta[1]     1   0.75
## beta[2]     1   0.62
## beta[3]     1   0.13
## beta[4]     1   0.36
## gamma       1   0.75
## theta       1   0.52
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    5110    7449    9000    9000    5453    3989    1248

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.97    1.36    0.90   -0.56   -0.72   -1.74    1.22 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.72
## beta[1]     1   0.27
## beta[2]     1   0.25
## beta[3]     1   0.38
## beta[4]     1   0.42
## gamma       1   0.16
## theta       1   0.20
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    4892    6488    8290    8695    5131    4037    2767

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.56   -0.60    0.56    0.77    0.23    0.56    0.58 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.18
## beta[1]     1   0.45
## beta[2]     1   0.76
## beta[3]     1   0.09
## beta[4]     1   0.28
## gamma       1   0.37
## theta       1   0.16
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    5138    6517    8150    9000    5225    4202    1330

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -1.55    0.48    1.56   -0.27    0.28    1.58   -1.79 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.58
## beta[1]     1   0.61
## beta[2]     1   0.24
## beta[3]     1   0.83
## beta[4]     1   0.85
## gamma       1   0.97
## theta       1   0.64
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    5213    7603    9000    9000    5615    3104     472

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -2.02   -1.62   -1.44   -1.73    0.22    3.80   -3.35 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.16
## beta[1]     1   0.59
## beta[2]     1   0.62
## beta[3]     1   0.30
## beta[4]     1   0.95
## gamma     901   0.38
## theta       1   0.60
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    4224    7568    9000    9000    5564    2169     204

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.53   -2.22    0.87   -0.70    0.18    1.16   -0.70 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.60
## beta[1]     1   0.22
## beta[2]     1   0.25
## beta[3]     1   0.49
## beta[4]     1   0.44
## gamma       1   0.83
## theta       1   0.54
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    5107    7594    9000    9000    5392    4375    1319

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.28    0.45   -1.01    0.56    0.70   -0.19   -0.01 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.37
## beta[1]     1   0.26
## beta[2]     1   0.80
## beta[3]     1   0.83
## beta[4]     1   0.53
## gamma       1   0.45
## theta       1   0.72
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    5024    7968    8951    9794    5237    4396     978

RMW Inverse Gamma mixing

ParamTableIGam <- expand.grid("Hyp1Gam" = c(4, 1),
                             "PriorMeanCV" = c(1.5, 5),
                             "PriorCV" = c("TruncExp", "Pareto"))
ParamTableIGam$Hyp2Gam = with(ParamTableIGam, 
                             ifelse(Hyp1Gam == 4, 1, ifelse(Hyp1Gam == 1, 1, 0.1)))
set.seed(103)
ParamTableIGam$seed = sample(1:1e6, size = nrow(ParamTableIGam))
cat("Parameter values \n")
## Parameter values
print(ParamTableIGam); cat("\n")
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed
## 1       4         1.5 TruncExp       1 215942
## 2       1         1.5 TruncExp       1  63169
## 3       4         5.0 TruncExp       1 521826
## 4       1         5.0 TruncExp       1 503616
## 5       4         1.5   Pareto       1 120486
## 6       1         1.5   Pareto       1  87275
## 7       4         5.0   Pareto       1 433560
## 8       1         5.0   Pareto       1 193080
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainIGam <- parLapply(cl, as.list(1:nrow(ParamTableIGam)), 
                      fun = RunRMWreg_MCMC,
                      ParamTable = ParamTableIGam, 
                      Mixing = "InvGamma", BaseModel = "Weibull",
                      N = N * 2, Thin = Thin * 2, Burn = Burn * 2, 
                      Time = Time, Event = Event, DesignMat = DesignMat, 
                      PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
                      StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)
ChainIGam_mcmc = lapply(ChainIGam, function(x) {mcmc(cbind(x$beta, x$gam, x$theta))})
ChainIGam_mcmc = lapply(ChainIGam_mcmc, function(x)
                                      {
                                        colnames(x) <- titles
                                        return(x)
                                      })
Plot = lapply(ChainIGam_mcmc, PlotDiag, titles = titles)

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    1.16   -1.21   -1.46   -0.55    0.52    0.93   -1.20 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.27
## beta[1]     1   0.48
## beta[2]     1   0.21
## beta[3]     1   0.88
## beta[4]     1   0.60
## gamma       1   0.35
## theta       1   0.45
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     315    9000    7819    8420    4446     456     114

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.36   -0.79   -0.60   -0.46   -0.69    0.28   -0.29 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.91
## beta[1]     1   0.42
## beta[2]     1   0.19
## beta[3]     1   0.86
## beta[4]     1   0.32
## gamma       1   0.87
## theta       1   0.93
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     331    8629    9000    9000    4967     529     187

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -1.51    2.10    1.39   -0.13    0.58   -1.68    1.11 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]  1801   0.56
## beta[1]   901   0.19
## beta[2]     1   0.34
## beta[3]     1   0.43
## beta[4]     1   0.50
## gamma    1801   0.18
## theta    1801   0.73
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     360    9000    9000    9000    4849     427     193

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.93   -1.70   -1.02    1.53    0.91    1.05   -1.19 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.28
## beta[1]     1   0.23
## beta[2]     1   0.66
## beta[3]     1   0.19
## beta[4]     1   0.68
## gamma       1   0.44
## theta       1   0.20
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     241    8055    9000    8020    4962     461     146

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.74    0.13   -0.55   -1.01   -1.39   -0.80    0.77 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.25
## beta[1]     1   0.61
## beta[2]     1   0.60
## beta[3]     1   0.39
## beta[4]     1   0.61
## gamma       1   0.36
## theta       1   0.09
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     230    9000    9000    9000    4620     474     111

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -1.20    1.30    1.38   -0.35    0.21   -1.45    0.81 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.85
## beta[1]     1   0.49
## beta[2]     1   0.44
## beta[3]     1   0.66
## beta[4]     1   0.20
## gamma       1   0.48
## theta       1   0.86
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     104    5578    9000    9000    5000     267      63

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    1.64    0.44    0.49    0.01   -0.72    1.06   -1.16 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]    NA   0.00
## beta[1]    NA   0.00
## beta[2]     1   0.83
## beta[3]     1   0.84
## beta[4]     1   0.21
## gamma      NA   0.00
## theta      NA   0.00
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##      42    5087    9000    8630    4750     251      12

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.96   -1.37   -1.31   -0.51   -0.43    0.94   -1.08 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]  1801   0.08
## beta[1]     1   0.52
## beta[2]     1   0.29
## beta[3]     1   0.74
## beta[4]     1   0.92
## gamma    1801   0.06
## theta       1   0.08
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     295    9000    9000    8105    5115     553     181

RMW Inverse Gaussian mixing

ParamTableIGauss <- expand.grid("Hyp1Gam" = c(4, 1),
                             "PriorMeanCV" = c(1.5, 5),
                             "PriorCV" = c("TruncExp", "Pareto"))
ParamTableIGauss$Hyp2Gam = with(ParamTableIGauss, 
                             ifelse(Hyp1Gam == 4, 1, ifelse(Hyp1Gam == 1, 1, 0.1)))
set.seed(104)
ParamTableIGauss$seed = sample(1:1e6, size = nrow(ParamTableIGauss))
cat("Parameter values \n")
## Parameter values
print(ParamTableIGauss); cat("\n")
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed
## 1       4         1.5 TruncExp       1 364452
## 2       1         1.5 TruncExp       1 772498
## 3       4         5.0 TruncExp       1 734878
## 4       1         5.0 TruncExp       1 972528
## 5       4         1.5   Pareto       1 740140
## 6       1         1.5   Pareto       1 200713
## 7       4         5.0   Pareto       1 377307
## 8       1         5.0   Pareto       1 247017
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainIGauss <- parLapply(cl, as.list(1:nrow(ParamTableIGauss)), 
                      fun = RunRMWreg_MCMC,
                      ParamTable = ParamTableIGauss, 
                      Mixing = "InvGauss", BaseModel = "Weibull",
                      N = 2*N, Thin = 2*Thin, Burn = 2*Burn, 
                      Time = Time, Event = Event, DesignMat = DesignMat, 
                      PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
                      StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)
ChainIGauss_mcmc = lapply(ChainIGauss, function(x) {mcmc(cbind(x$beta, x$gam, x$theta))})
ChainIGauss_mcmc = lapply(ChainIGauss_mcmc, function(x)
                                      {
                                        colnames(x) <- titles
                                        return(x)
                                      })
Plot = lapply(ChainIGauss_mcmc, PlotDiag, titles = titles)

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.83    0.79    0.80   -0.77   -0.96    0.11    0.61 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.24
## beta[1]     1   0.07
## beta[2]     1   0.88
## beta[3]     1   0.95
## beta[4]     1   0.53
## gamma       1   0.19
## theta       1   0.06
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     298    5342    7007    7754    1096     217     134

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.70    0.25    0.52    1.06    0.36   -1.16   -1.60 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]    NA   0.00
## beta[1]     1   0.09
## beta[2]     1   0.12
## beta[3]     1   0.66
## beta[4]     1   0.12
## gamma       1   0.06
## theta       1   0.15
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     296    3266    7750    8240    1040     142     115

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.03    0.38   -0.91    1.63   -0.18    0.26    0.05 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.05
## beta[1]     1   0.14
## beta[2]     1   0.91
## beta[3]     1   0.14
## beta[4]     1   0.79
## gamma       1   0.12
## theta    1801   0.07
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     843    6211    8506    9000    4218     531     229

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -1.06    1.23    0.37   -1.27   -0.53   -0.99   -0.69 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.27
## beta[1]     1   0.33
## beta[2]     1   0.61
## beta[3]     1   0.44
## beta[4]     1   0.51
## gamma       1   0.49
## theta       1   0.72
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     597    6137    8700    9000    3971     442     341

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    1.49   -0.86   -0.31    1.55   -0.78    0.42   -0.36 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.37
## beta[1]     1   0.36
## beta[2]     1   0.07
## beta[3]     1   0.46
## beta[4]     1   0.64
## gamma       1   0.49
## theta       1   0.69
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     478    5766    9000    9000    4171     469     424

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.51    0.32   -0.71   -0.13    0.31   -0.41   -0.26 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.90
## beta[1]     1   0.57
## beta[2]     1   0.20
## beta[3]     1   0.94
## beta[4]     1   0.52
## gamma       1   0.82
## theta       1   0.80
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     323    4792    9000    8370    4278     372     275

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.50   -0.12   -0.29    0.96   -1.17    0.14   -0.61 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.63
## beta[1]     1   0.48
## beta[2]     1   0.34
## beta[3]     1   0.32
## beta[4]     1   0.22
## gamma       1   0.49
## theta       1   0.63
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     559    6760    8362    9000    4113     447     360

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    2.12   -1.77    0.65   -0.23    0.81    2.17    1.33 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.26
## beta[1]     1   0.63
## beta[2]     1   0.40
## beta[3]     1   0.85
## beta[4]     1   0.60
## gamma     901   0.09
## theta    2701   0.44
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##     372    5317    8534    9000    4180     310     194

RMW Log Normal mixing

ParamTableLN <- expand.grid("Hyp1Gam" = c(4, 1),
                             "PriorMeanCV" = c(1.5, 5),
                             "PriorCV" = c("TruncExp", "Pareto"))
ParamTableLN$Hyp2Gam = with(ParamTableLN, 
                             ifelse(Hyp1Gam == 4, 1, ifelse(Hyp1Gam == 1, 1, 0.1)))
cat("Parameter values \n")
## Parameter values
print(ParamTableLN); cat("\n")
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam
## 1       4         1.5 TruncExp       1
## 2       1         1.5 TruncExp       1
## 3       4         5.0 TruncExp       1
## 4       1         5.0 TruncExp       1
## 5       4         1.5   Pareto       1
## 6       1         1.5   Pareto       1
## 7       4         5.0   Pareto       1
## 8       1         5.0   Pareto       1
set.seed(105)
ParamTableLN$seed = sample(1:1e6, size = nrow(ParamTableLN))
cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ChainLN <- parLapply(cl, as.list(1:nrow(ParamTableLN)), 
                      fun = RunRMWreg_MCMC,
                      ParamTable = ParamTableLN, 
                      Mixing = "LogNormal", BaseModel = "Weibull",
                      N = 2*N, Thin = 2*Thin, Burn = 2*Burn, 
                      Time = Time, Event = Event, DesignMat = DesignMat, 
                      PrintProgress = FALSE, StoreAdapt = TRUE, lambdaPeriod = 1,
                      StoreChains = TRUE, StoreDir = StoreDir, DataName = "CP")
stopCluster(cl)
ChainLN_mcmc = lapply(ChainLN, function(x) {mcmc(cbind(x$beta, x$gam, x$theta))})
ChainLN_mcmc = lapply(ChainLN_mcmc, function(x)
                                      {
                                        colnames(x) <- titles
                                        return(x)
                                      })
Plot = lapply(ChainLN_mcmc, PlotDiag, titles = titles)

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.40   -1.89   -0.99   -0.84   -0.92    1.88    1.02 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.68
## beta[1]     1   0.07
## beta[2]     1   0.58
## beta[3]     1   0.07
## beta[4]     1   0.64
## gamma       1   0.13
## theta       1   0.24
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    3359    5908    8689    9000    3504    1044     980

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.16    0.38    0.13   -1.55    0.11   -0.36    0.01 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.82
## beta[1]     1   0.28
## beta[2]     1   0.25
## beta[3]     1   0.29
## beta[4]     1   0.88
## gamma       1   0.18
## theta       1   0.29
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    3418    4761    8380    9000    3581     841     750

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.87   -0.83    1.16    1.34    0.67    0.02   -0.16 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.37
## beta[1]     1   0.33
## beta[2]     1   0.79
## beta[3]     1   0.43
## beta[4]     1   0.36
## gamma       1   0.94
## theta       1   0.89
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    3157    6741    8267    9000    3302    1055     978

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    1.54    0.16   -0.71    0.82   -1.65    0.53    0.68 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.85
## beta[1]     1   0.13
## beta[2]     1   0.40
## beta[3]     1   0.19
## beta[4]     1   0.48
## gamma       1   0.49
## theta       1   0.69
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    2918    5410    8573    9000    3331     844     856

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -1.10    1.68   -0.59   -0.87    1.85   -2.09   -2.14 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.27
## beta[1]     1   0.50
## beta[2]     1   0.61
## beta[3]     1   0.70
## beta[4]     1   0.11
## gamma       1   0.13
## theta       1   0.07
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    3776    6280    8587    9000    3848     851     846

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.34    0.11   -0.18    1.17   -0.22    0.53    0.51 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.81
## beta[1]     1   0.52
## beta[2]     1   0.47
## beta[3]     1   0.39
## beta[4]     1   0.90
## gamma       1   0.73
## theta       1   0.76
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    4043    3612    8518    8605    3915     581     455

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    0.90    2.24    0.50    0.33   -0.26   -1.54   -0.84 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.32
## beta[1]  1801   0.21
## beta[2]     1   0.12
## beta[3]     1   0.22
## beta[4]     1   0.49
## gamma       1   0.08
## theta       1   0.34
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    3177    6390    8360    8768    3550     871     730

## Geweke diagnostics Z scores 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##   -0.52    1.77    1.76    0.59    0.86   -0.65   -0.79 
## 
## Heidelberger and Welch's convergence diagnostic 
##         start pvalue
## beta[0]     1   0.13
## beta[1]     1   0.28
## beta[2]     1   0.40
## beta[3]     1   0.56
## beta[4]     1   0.13
## gamma       1   0.14
## theta       1   0.15
## 
## Effective Sample Size 
## beta[0] beta[1] beta[2] beta[3] beta[4]   gamma   theta 
##    3544    5673    8501    9000    3650     724     730

Table 5

cv<-function(gam,theta,mixing,ratio=FALSE)
{
    if(mixing=="Gamma")
    {
        cv.star2=exp(lgamma(theta)+lgamma(theta-2/gam)-2*lgamma(theta-1/gam))-1
        dcv.star2=exp(lgamma(theta)+lgamma(theta-2/gam)-2*lgamma(theta-1/gam))*(digamma(theta)+digamma(theta-2/gam)-2*digamma(theta-1/gam))
    }
    if(mixing=="InvGamma")
    {
        cv.star2=exp(lgamma(theta)+lgamma(theta+2/gam)-2*lgamma(theta+1/gam))-1
        dcv.star2=exp(lgamma(theta)+lgamma(theta+2/gam)-2*lgamma(theta+1/gam))*(digamma(theta)+digamma(theta+2/gam)-2*digamma(theta+1/gam))
    }
    if(mixing=="InvGauss")
    {
        K1=besselK(x=1/theta,nu=-(1/gam+1/2)); K2=besselK(x=1/theta,nu=-(2/gam+1/2))
        
        #Using the asymptotic expansion as theta->0
        cv.star2.aux=sqrt(pi/2)*(exp(-1/theta)*theta^(1/2))*(K2/K1)*(1/K1)-1
        aux1=sqrt(pi/2)*theta^(-3/2)*exp(-1/theta)/K1
        aux2=(K2/K1)+besselK(x=1/theta,nu=-(2/gam-1/2))/K1-2*(K2/K1)*(1/K1)*besselK(x=1/theta,nu=-(1/gam-1/2))
        dcv.star2.aux=aux1*aux2

        cv.star2=ifelse(theta<0.0015,0,cv.star2.aux)
        dcv.star2=ifelse(theta<0.0015,0,dcv.star2.aux)

    }
    if(mixing=="LogNormal")
    {
        cv.star2=exp(theta)-1
        dcv.star2=exp(-2*log(gam)+theta)
    }

    aux0=(lgamma(1+2/gam)-2*lgamma(1+1/gam)); aux1=exp(aux0)-1
    cv=sqrt(exp(aux0+log(cv.star2))+aux1)

    if(ratio==FALSE){return(cv)}
    if(ratio==TRUE){return(cv/sqrt(aux1))}
    if(ratio=="WEI"){return(sqrt(aux1))}
}

cv.Weibull<-function(gam)
{
    aux=sqrt(exp(lgamma(1+2/gam)-2*lgamma(1+1/gam))-1)
    return(aux)
}

RcvGam = matrix(0, ncol = 3, nrow = 8)
for(i in 1:nrow(ParamTableGam))
{
  aux = mcmc(cv(ChainGam[[i]]$gam,ChainGam[[i]]$theta,mixing="Gamma",ratio=TRUE))
  RcvGam[i, 1] = median(aux)
  RcvGam[i, 2:3] = HPDinterval(aux)
}
colnames(RcvGam) <- c("Ecv", "HPDEcvlower", "HPDEcvupper")
RcvIGam = matrix(0, ncol = 3, nrow = 8)
for(i in 1:nrow(ParamTableIGam))
{
  aux = mcmc(cv(ChainIGam[[i]]$gam,ChainIGam[[i]]$theta,mixing="InvGamma",ratio=TRUE))
  RcvIGam[i, 1] = median(aux)
  RcvIGam[i, 2:3] = HPDinterval(aux)
}
colnames(RcvIGam) <- c("Ecv", "HPDEcvlower", "HPDEcvupper")
RcvIGauss = matrix(0, ncol = 3, nrow = 8)
for(i in 1:nrow(ParamTableIGauss))
{
  aux = mcmc(cv(ChainIGauss[[i]]$gam,ChainIGauss[[i]]$theta,mixing="InvGauss",ratio=TRUE))
  RcvIGauss[i, 1] = median(aux)
  RcvIGauss[i, 2:3] = HPDinterval(aux)
}
colnames(RcvIGauss) <- c("Ecv", "HPDEcvlower", "HPDEcvupper")
RcvLN = matrix(0, ncol = 3, nrow = 8)
for(i in 1:nrow(ParamTableLN))
{
  aux = mcmc(cv(ChainLN[[i]]$gam,ChainLN[[i]]$theta,mixing="LogNormal",ratio=TRUE))
  RcvLN[i, 1] = median(aux)
  RcvLN[i, 2:3] = HPDinterval(aux)
}
colnames(RcvLN) <- c("Ecv", "HPDEcvlower", "HPDEcvupper")

cat("Estimated ratios for each model \n")
## Estimated ratios for each model
cbind(ParamTableGam, RcvGam)
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed      Ecv HPDEcvlower
## 1       4         1.5 TruncExp       1 571629 2.391995    1.226446
## 2       1         1.5 TruncExp       1 492883 2.322744    1.190385
## 3       4         5.0 TruncExp       1 783694 6.986370    1.589244
## 4       1         5.0 TruncExp       1 540340 6.943675    1.448344
## 5       4         1.5   Pareto       1  88002 2.329099    1.087620
## 6       1         1.5   Pareto       1 420628 2.240850    1.040911
## 7       4         5.0   Pareto       1 976328 4.344738    1.146542
## 8       1         5.0   Pareto       1 330313 4.070118    1.159338
##   HPDEcvupper
## 1    4.422618
## 2    4.235218
## 3   19.874946
## 4   19.789993
## 5    6.134094
## 6    5.767696
## 7   29.566035
## 8   27.916060
cbind(ParamTableIGam, RcvIGam)
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed      Ecv HPDEcvlower
## 1       4         1.5 TruncExp       1 215942 1.422416    1.243369
## 2       1         1.5 TruncExp       1  63169 1.407057    1.209227
## 3       4         5.0 TruncExp       1 521826 1.430164    1.237122
## 4       1         5.0 TruncExp       1 503616 1.403782    1.194593
## 5       4         1.5   Pareto       1 120486 1.399151    1.192560
## 6       1         1.5   Pareto       1  87275 1.364267    1.131292
## 7       4         5.0   Pareto       1 433560 1.422207    1.198182
## 8       1         5.0   Pareto       1 193080 1.400980    1.201202
##   HPDEcvupper
## 1    1.554756
## 2    1.553597
## 3    1.556096
## 4    1.551413
## 5    1.551991
## 6    1.539696
## 7    1.570126
## 8    1.546437
cbind(ParamTableIGauss, RcvIGauss)
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed      Ecv HPDEcvlower
## 1       4         1.5 TruncExp       1 364452 1.671872    1.467143
## 2       1         1.5 TruncExp       1 772498 1.662514    1.452708
## 3       4         5.0 TruncExp       1 734878 1.676073    1.471130
## 4       1         5.0 TruncExp       1 972528 1.658172    1.435204
## 5       4         1.5   Pareto       1 740140 1.646229    1.395405
## 6       1         1.5   Pareto       1 200713 1.620578    1.352898
## 7       4         5.0   Pareto       1 377307 1.669073    1.436302
## 8       1         5.0   Pareto       1 247017 1.637123    1.390757
##   HPDEcvupper
## 1    1.830358
## 2    1.831276
## 3    1.848358
## 4    1.838652
## 5    1.831929
## 6    1.817909
## 7    1.837317
## 8    1.829701
cbind(ParamTableLN, RcvLN)
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed      Ecv HPDEcvlower
## 1       4         1.5 TruncExp       1  97915 2.384842    1.533376
## 2       1         1.5 TruncExp       1 995633 2.238037    1.457315
## 3       4         5.0 TruncExp       1 320552 2.559325    1.689047
## 4       1         5.0 TruncExp       1 397726 2.431954    1.611196
## 5       4         1.5   Pareto       1 553224 2.283492    1.519334
## 6       1         1.5   Pareto       1 795749 2.075842    1.218768
## 7       4         5.0   Pareto       1 902245 2.434995    1.578776
## 8       1         5.0   Pareto       1 827927 2.247022    1.415542
##   HPDEcvupper
## 1    3.169264
## 2    3.126312
## 3    3.388848
## 4    3.289663
## 5    3.159489
## 6    2.969993
## 7    3.291863
## 8    3.152258
cat("Estimated true ratio w.r.t. Weibull model (Gamma model only with gam~Gamma(1,4))")
## Estimated true ratio w.r.t. Weibull model (Gamma model only with gam~Gamma(1,4))
aux = cv(ChainGam[[1]]$gam,ChainGam[[1]]$theta,mixing="Gamma") / cv.Weibull(ChainNone[[1]]$gam)
round(cbind(median(aux), HPDinterval(mcmc(aux))),2)
##         lower upper
## V1 2.08  1.08  3.74

Model comparison

Pseudo ML

CDnone = lapply(ChainNone, FUN = RMWreg_CaseDeletion, 
               Time = Time, Event = Event, DesignMat = DesignMat, 
               Mixing = "None", BaseModel = "Weibull")
ParamTableNone$logPsML = sapply(CDnone, FUN = function(x) sum(x[,1]))

CDexp = lapply(ChainExp, FUN = RMWreg_CaseDeletion, 
               Time = Time, Event = Event, DesignMat = DesignMat, 
               Mixing = "Exponential", BaseModel = "Weibull")
ParamTableExp$logPsML = sapply(CDexp, FUN = function(x) sum(x[,1]))

CDgam = lapply(ChainGam, FUN = RMWreg_CaseDeletion, 
               Time = Time, Event = Event, DesignMat = DesignMat, 
               Mixing = "Gamma", BaseModel = "Weibull")
ParamTableGam$logPsML = sapply(CDgam, FUN = function(x) sum(x[,1]))

CDigam = lapply(ChainIGam, FUN = RMWreg_CaseDeletion, 
                Time = Time, Event = Event, DesignMat = DesignMat, 
                Mixing = "InvGamma", BaseModel = "Weibull")
ParamTableIGam$logPsML = sapply(CDigam, FUN = function(x) sum(x[,1]))

CDigauss = lapply(ChainIGauss, FUN = RMWreg_CaseDeletion, 
                  Time = Time, Event = Event, DesignMat = DesignMat, 
                  Mixing = "InvGauss", BaseModel = "Weibull")
ParamTableIGauss$logPsML = sapply(CDigauss, FUN = function(x) sum(x[,1]))

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
CDln = parLapply(cl, ChainLN, fun = RMWreg_CaseDeletion, 
                 Time = Time, Event = Event, DesignMat = DesignMat, 
                 Mixing = "LogNormal", BaseModel = "Weibull")
stopCluster(cl)
ParamTableLN$logPsML = sapply(CDln, FUN = function(x) sum(x[,1]))

ML

RunRMWreg_logML <- function(Row, ParamTable, Chains,
                            Mixing, BaseModel,
                            Time, Event, DesignMat, 
                            N, Thin, Burn, lambdaPeriod = 1)
{
  require(RMWreg)
  
  if(Mixing %in% c("None", "Exponential"))
  {
    logML = RMWreg_logML(Chains[[Row]], Time, Event, DesignMat, 
                         Hyp1Gam = ParamTable$Hyp1Gam[Row], Hyp2Gam = ParamTable$Hyp2Gam[Row],
                         Mixing = Mixing, BaseModel = BaseModel, 
                         N=N, Thin = Thin, Burn = Burn,
                         lambdaPeriod = lambdaPeriod)        
  }
  else
  {
    logML = RMWreg_logML(Chains[[Row]], Time, Event, DesignMat, 
                         PriorCV = as.character(ParamTable$PriorCV[Row]),
                         PriorMeanCV = ParamTable$PriorMeanCV[Row],
                         Hyp1Gam = ParamTable$Hyp1Gam[Row], Hyp2Gam = ParamTable$Hyp2Gam[Row],
                         Mixing = Mixing, BaseModel = BaseModel, 
                         N=N, Thin = Thin, Burn = Burn,
                         lambdaPeriod = lambdaPeriod)    
  }
  
  write.table(logML, paste0("logML_",Mixing,"_",Row,".txt"), row.names = FALSE, col.names = FALSE)

  return(logML)
  
}

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableNone$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableNone)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableNone, 
                                      Chains = ChainNone,
                                      Mixing = "None", 
                                      BaseModel = "Weibull",
                                      Time = Time, Event = Event,
                                      DesignMat = DesignMat, 
                                      N=N, Thin = Thin, Burn = Burn,
                                      lambdaPeriod = 1))
stopCluster(cl)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableExp$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableExp)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableExp, 
                                      Chains = ChainExp,
                                      Mixing = "Exponential",
                                      BaseModel = "Weibull",
                                      Time = Time, Event = Event,
                                      DesignMat = DesignMat, 
                                      N=N, Thin = Thin, Burn = Burn,
                                      lambdaPeriod = 1))
stopCluster(cl)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableGam$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableGam)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableGam, 
                                      Chains = ChainGam,
                                      Mixing = "Gamma", 
                                      BaseModel = "Weibull",
                                      Time = Time, Event = Event,
                                      DesignMat = DesignMat, 
                                      N=N, Thin = Thin, Burn = Burn,
                                      lambdaPeriod = 1))
stopCluster(cl)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableIGam$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableIGam)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableIGam, 
                                      Chains = ChainIGam,
                                      Mixing = "InvGamma", BaseModel = "Weibull",
                                      Time = Time, Event = Event, 
                                      DesignMat = DesignMat, 
                                      N=N, Thin = Thin, Burn = Burn,
                                      lambdaPeriod = 1))
stopCluster(cl)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableIGauss$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableIGauss)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableIGauss, 
                                      Chains = ChainIGauss,
                                      Mixing = "InvGauss", BaseModel = "Weibull",
                                      Time = Time, Event = Event, DesignMat = DesignMat, 
                                      N=N, Thin = Thin, Burn = Burn,
                                      lambdaPeriod = 1))
stopCluster(cl)

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableLN$logML <- simplify2array(parLapply(cl, as.list(1:nrow(ParamTableLN)), 
                                      fun = RunRMWreg_logML,
                                      ParamTable = ParamTableLN, 
                                      Chains = ChainLN,
                                      Mixing = "LogNormal", BaseModel = "Weibull",
                                      Time = Time, Event = Event, DesignMat = DesignMat, 
                                      N=N, Thin = Thin, Burn = Burn,
                                      lambdaPeriod = 1))
stopCluster(cl)

DIC

ParamTableNone$DIC <- sapply(ChainNone, FUN = RMWreg_DIC, 
                           Time = Time, Event = Event, DesignMat = DesignMat, 
                           Mixing = "None", BaseModel = "Weibull")

ParamTableExp$DIC <- sapply(ChainExp, FUN = RMWreg_DIC, 
                           Time = Time, Event = Event, DesignMat = DesignMat, 
                           Mixing = "Exponential", BaseModel = "Weibull")

ParamTableGam$DIC = sapply(ChainGam, FUN = RMWreg_DIC, 
                           Time = Time, Event = Event, DesignMat = DesignMat, 
                           Mixing = "Gamma", BaseModel = "Weibull")

ParamTableIGam$DIC = sapply(ChainIGam, FUN = RMWreg_DIC, 
                            Time = Time, Event = Event, DesignMat = DesignMat, 
                            Mixing = "InvGamma", BaseModel = "Weibull")

ParamTableIGauss$DIC = sapply(ChainIGauss, FUN = RMWreg_DIC, 
                              Time = Time, Event = Event, DesignMat = DesignMat, 
                              Mixing = "InvGauss", BaseModel = "Weibull")

cl <- makeCluster(getOption("cl.cores", detectCores()-1))
ParamTableLN$DIC = simplify2array(parLapply(cl, ChainLN, fun = RMWreg_DIC, 
                                            Time = Time, Event = Event, 
                                            DesignMat = DesignMat, 
                                            Mixing = "LogNormal", 
                                            BaseModel = "Weibull"))
stopCluster(cl)

Summary

cat("\n RMW no mixing \n")
## 
##  RMW no mixing
ParamTableNone
##   Hyp1Gam Hyp2Gam   seed   logPsML     logML
## 1       4       1 307767 -1235.611 -1240.090
## 2       1       1 257673 -1235.667 -1238.788
cat("\n RMW Exponential mixing \n")
## 
##  RMW Exponential mixing
ParamTableExp
##   Hyp1Gam Hyp2Gam   seed   logPsML     logML
## 1       4       1 372199 -1227.396 -1231.353
## 2       1       1  43825 -1227.506 -1230.724
cat("\n RMW Gamma mixing \n")
## 
##  RMW Gamma mixing
ParamTableGam
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed   logPsML     logML
## 1       4         1.5 TruncExp       1 571629 -1230.284 -1234.942
## 2       1         1.5 TruncExp       1 492883 -1230.294 -1233.972
## 3       4         5.0 TruncExp       1 783694 -1229.227 -1233.306
## 4       1         5.0 TruncExp       1 540340 -1229.285 -1232.612
## 5       4         1.5   Pareto       1  88002 -1230.321 -1235.412
## 6       1         1.5   Pareto       1 420628 -1230.487 -1234.411
## 7       4         5.0   Pareto       1 976328 -1229.537 -1234.149
## 8       1         5.0   Pareto       1 330313 -1229.655 -1233.316
cat("\n RMW Inverse Gamma mixing \n")
## 
##  RMW Inverse Gamma mixing
ParamTableIGam
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed   logPsML     logML
## 1       4         1.5 TruncExp       1 215942 -1230.997 -1231.651
## 2       1         1.5 TruncExp       1  63169 -1231.057 -1233.952
## 3       4         5.0 TruncExp       1 521826 -1230.985 -1235.228
## 4       1         5.0 TruncExp       1 503616 -1231.166 -1234.474
## 5       4         1.5   Pareto       1 120486 -1231.048 -1235.760
## 6       1         1.5   Pareto       1  87275 -1231.385 -1234.472
## 7       4         5.0   Pareto       1 433560 -1231.030 -1233.607
## 8       1         5.0   Pareto       1 193080 -1231.065 -1234.788
cat("\n RMW Inverse Gaussian mixing \n")
## 
##  RMW Inverse Gaussian mixing
ParamTableIGauss
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed   logPsML     logML
## 1       4         1.5 TruncExp       1 364452 -1228.718 -1233.280
## 2       1         1.5 TruncExp       1 772498 -1228.764 -1232.417
## 3       4         5.0 TruncExp       1 734878 -1228.781 -1233.141
## 4       1         5.0 TruncExp       1 972528 -1228.803 -1232.628
## 5       4         1.5   Pareto       1 740140 -1228.960 -1233.690
## 6       1         1.5   Pareto       1 200713 -1229.028 -1233.471
## 7       4         5.0   Pareto       1 377307 -1228.847 -1233.492
## 8       1         5.0   Pareto       1 247017 -1228.974 -1232.903
cat("\n RMW Log-Normal mixing \n")
## 
##  RMW Log-Normal mixing
ParamTableLN
##   Hyp1Gam PriorMeanCV  PriorCV Hyp2Gam   seed   logPsML     logML
## 1       4         1.5 TruncExp       1  97915 -1227.920 -1232.143
## 2       1         1.5 TruncExp       1 995633 -1228.143 -1231.776
## 3       4         5.0 TruncExp       1 320552 -1227.818 -1232.327
## 4       1         5.0 TruncExp       1 397726 -1228.007 -1232.309
## 5       4         1.5   Pareto       1 553224 -1228.010 -1232.904
## 6       1         1.5   Pareto       1 795749 -1228.567 -1232.309
## 7       4         5.0   Pareto       1 902245 -1227.877 -1232.262
## 8       1         5.0   Pareto       1 827927 -1228.286 -1231.900

Figure 9

PosteriorSummary <- function(Column, ChainNone, ChainExp, ChainGam,
                             ChainIGam, ChainIGauss, ChainLN)
{
  Median = c(unlist(lapply(ChainNone, function(x) median(x[,Column]))),
             unlist(lapply(ChainExp, function(x) median(x[,Column]))),
             unlist(lapply(ChainGam, function(x) median(x[,Column]))),
             unlist(lapply(ChainIGam, function(x) median(x[,Column]))),
             unlist(lapply(ChainIGauss, function(x) median(x[,Column]))),
             unlist(lapply(ChainLN, function(x) median(x[,Column]))))  
  HPD = rbind(t(sapply(ChainNone, function(x) HPDinterval(x[,Column]))),
              t(sapply(ChainExp, function(x) HPDinterval(x[,Column]))),
              t(sapply(ChainGam, function(x) HPDinterval(x[,Column]))),
              t(sapply(ChainIGam, function(x) HPDinterval(x[,Column]))),
              t(sapply(ChainIGauss, function(x) HPDinterval(x[,Column]))),
              t(sapply(ChainLN, function(x) HPDinterval(x[,Column]))))
  
  return(cbind(Median, HPD))
}

plotHPD <- function(hpd, ...)
{
    plotCI(x = hpd[,1], uiw = hpd[,3]-hpd[,1], liw = hpd[,1]-hpd[,2],
           xaxt = "n", lwd = 2, xlab = "Mixing-Prior", gap = 1.5, ...)
}

param.name = c(expression(beta[0]), expression(beta[1]), expression(beta[2]), expression(beta[3]),
               expression(beta[4]), expression(gamma), expression(theta))
titles = c(expression("Intercept"), expression("No impairments"), 
           expression("1 impairment"), expression("2 impairments"), expression("Birth weight"),"")

par(mfrow = c(6,1))
par(mar = c(7.5,11,6,1), mgp = c(6, 2, 0))
for(Column in 1:6)
{
  PS = PosteriorSummary(Column, ChainNone_mcmc, ChainExp_mcmc, ChainGam_mcmc,
                                ChainIGam_mcmc, ChainIGauss_mcmc, ChainLN_mcmc)  
  if(Column == 1)  
  { 
    plotHPD(PS, main = titles[Column], ylab = param.name[Column],
            cex.lab = 2.5, cex.axis = 2.5, cex.main = 4, ylim = c(2.5,4.8))
  }
  else
  {
    plotHPD(PS, main = titles[Column], ylab = param.name[Column],
            cex.lab = 2.5, cex.axis = 2.5, cex.main = 4)    
  }
  
  abline(v = c(2.5,4.5,8.5,12.5,16.5,20.5,24.5,28.5,32.5), lty = 1, lwd = 4)
  abline(v = c(6.5,10.5,14.5,18.5,22.5,26.5,30.5,34.5), lty = 3)
  axis(side = 1, at = c(1.5,3.5,6.5,10.5,14.5,18.5,22.5,26.5,30.5,34.5),
       labels = c("None", "Exp", "Gamma-TExp", "Gam-Pare", "IGam-TExp", "IGam-Pareto", 
                  "IGauss-TEx", "IGauss-Pare", "LN-TExp", "LN-Pareto"), cex.axis = 2)
  if(Column == 1)
  {
    text(seq(5.5,36,by=2), rep(4.5,16),
         rep(c(expression(paste("E(",c[v],")=1.5")), 
               expression(paste("E(",c[v],")=5"))),8), cex = 1.5)
  }
}
## Note: no visible binding for global variable 'v' 
## Note: no visible binding for global variable 'v'

Figure 10

PlotModelComparison <- function(Row, ParamValue, Mains,
                         ParamTableNone, ParamTableExp, ParamTableGam, 
                         ParamTableIGam, ParamTableIGauss, ParamTableLN, ...)
{
  logBF = c(ParamTableNone$logML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]], 
             ParamTableExp$logML[ParamTableExp$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
             ParamTableGam$logML[ParamTableGam$PriorCV == "TruncExp" &
                                 ParamTableGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                 ParamTableGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
             ParamTableIGam$logML[ParamTableIGam$PriorCV == "TruncExp" & 
                                  ParamTableIGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                  ParamTableIGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
             ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "TruncExp" & 
                                    ParamTableIGauss$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                    ParamTableIGauss$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
             ParamTableLN$logML[ParamTableLN$PriorCV == "TruncExp" & 
                                ParamTableLN$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                ParamTableLN$Hyp1Gam == ParamValue$Hyp1Gam[Row]]) -
    ParamTableNone$logML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]]
  
  logBF_1 = c(ParamTableGam$logML[ParamTableGam$PriorCV == "Pareto" &
                                 ParamTableGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                 ParamTableGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
           ParamTableIGam$logML[ParamTableIGam$PriorCV == "Pareto" & 
                                  ParamTableIGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                  ParamTableIGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
           ParamTableIGauss$logML[ParamTableIGauss$PriorCV == "Pareto" & 
                                    ParamTableIGauss$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                    ParamTableIGauss$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
           ParamTableLN$logML[ParamTableLN$PriorCV == "Pareto" & 
                                ParamTableLN$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                ParamTableLN$Hyp1Gam == ParamValue$Hyp1Gam[Row]]) -
    ParamTableNone$logML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]]
  
    logPsBF = c(ParamTableNone$logPsML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]], 
             ParamTableExp$logPsML[ParamTableExp$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
             ParamTableGam$logPsML[ParamTableGam$PriorCV == "TruncExp" &
                                 ParamTableGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                 ParamTableGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
             ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "TruncExp" & 
                                  ParamTableIGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                  ParamTableIGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
             ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "TruncExp" & 
                                    ParamTableIGauss$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                    ParamTableIGauss$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
             ParamTableLN$logPsML[ParamTableLN$PriorCV == "TruncExp" & 
                                ParamTableLN$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                ParamTableLN$Hyp1Gam == ParamValue$Hyp1Gam[Row]]) -
    ParamTableNone$logPsML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]]
  
  logPsBF_1 = c(ParamTableGam$logPsML[ParamTableGam$PriorCV == "Pareto" &
                                 ParamTableGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                 ParamTableGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
           ParamTableIGam$logPsML[ParamTableIGam$PriorCV == "Pareto" & 
                                  ParamTableIGam$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                  ParamTableIGam$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
           ParamTableIGauss$logPsML[ParamTableIGauss$PriorCV == "Pareto" & 
                                    ParamTableIGauss$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                    ParamTableIGauss$Hyp1Gam == ParamValue$Hyp1Gam[Row]],
           ParamTableLN$logPsML[ParamTableLN$PriorCV == "Pareto" & 
                                ParamTableLN$PriorMeanCV == ParamValue$PriorMeanCV[Row] &
                                ParamTableLN$Hyp1Gam == ParamValue$Hyp1Gam[Row]]) -
    ParamTableNone$logPsML[ParamTableNone$Hyp1Gam == ParamValue$Hyp1Gam[Row]]
  
  plot(logBF, logPsBF, pch=c(4,8,0,5,1,2), cex=rep(3,6), bty = "n", 
       ylab = "Log-Pseudo Bayes Factors", xlab = "log-Bayes Factors", 
       ylim = c(0,9), xlim = c(0,9), main = Mains[Row], ...) 
  points(logBF_1, logPsBF_1, pch=c(15,18,16,17), cex=c(3,4.5,3,3), ...)
  
  if(Row == 4)
  {
    legend('bottomright', c("Weibull", "RMWEXP", "RMWGAM", "RMWIGAM", "RMWIGAUSS", "RMWLN"),
           pch = c(4,8,0,5,1,2), col = "black", cex = 1.5, bty = "n")
  }
}

ParamValue = expand.grid("Hyp1Gam" = c(4, 1),
                          "PriorMeanCV" = c(1.5, 5))
Mains = rep(c(expression(paste(gamma,"~Gamma(",4,",",1,")")), 
          expression(paste(gamma,"~Gamma(",1,",",1,")"))),2)

par(mar=c(5,5,4,0.5))
par(mfrow=c(2,2))
for(Row in 1:4)
{
  PlotModelComparison(Row, ParamValue, Mains,
                      ParamTableNone, ParamTableExp, ParamTableGam, 
                      ParamTableIGam, ParamTableIGauss, ParamTableLN,
                      cex.main = 3, cex.axis = 2, cex.lab = 2.5, lwd = 2)
}

Figure 11

# Outlier detection
BFexp = RMWreg_BFoutlier(ChainExp[[1]], RefLambda = Event + (1-Event)*1/2,
                         Time, Event, DesignMat,
                         Mixing = "Exponential", BaseModel = "Weibull")
par(mar = c(5, 5, 4, 2) + 0.1)
plot(BFexp, type = "l", ylim = c(0, 2), main = "",
     ylab = "Bayes Factors", xlab = "Patient", 
     cex.main = 3, cex.axis = 2, cex.lab = 2.5)
abline(h = 1, lty = 2)

Session information

sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.3
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] compiler  parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
## [1] KMsurv_0.1-5     data.table_1.9.6 gplots_3.0.1     coda_0.18-1     
## [5] RMWreg_0.0.18    knitr_1.15.1    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.8        magrittr_1.5       lattice_0.20-34   
##  [4] stringr_1.1.0      highr_0.6          caTools_1.17.1    
##  [7] tools_3.3.2        grid_3.3.2         KernSmooth_2.23-15
## [10] htmltools_0.3.5    gtools_3.5.0       yaml_2.1.14       
## [13] rprojroot_1.1      digest_0.6.10      bitops_1.0-6      
## [16] evaluate_0.10      rmarkdown_1.2      gdata_2.17.0      
## [19] stringi_1.1.2      backports_1.0.4    chron_2.3-47